#
# data <- readRDS('data_sim_p_1000.rds')
# for(x in names(data)) assign(x, data[[x]], envir = globalenv())
# U_full <- data$U
# rm(data)
#
# S.data.time <- S.data$obs
# S.data.time.log.sum <- sum(log(S.data.time))
# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, phi1, phi2, phi3) (phi1 )/(1+exp((phi2-t)/phi3))
#=======================================#
p <- 1000
parameter <- list(sigma2 = .05^2,
mu = c(0.9,90,5),
omega2 = c(0.005, 40, 1),
bara = 90,
barb = 30,
baralpha = 0.5,
beta = rep(0,p))
parameter$beta[1:4] <- c(1,2,-1,-2)
#=======================================#
t <- seq(60,120, length.out = 10) #time values
set.seed(123)
G <- 40 ; ng = 4
link = m
dt <- create_JLS_HD_data(G, ng, t, m, link, parameter)
var.true <- dt$var.true
a <- var.true$a ; var.true$a <- NULL #a fixé (et retiré des variables latentes)
S.data <- dt$survival
U_full <- dt$U
Y <- do.call(get_obs, var.true) + rnorm(n, 0, sqrt(parameter$sigma2))
S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model(
function(sigma2, ...) -n/(2*sigma2),
function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
# === Variable Latente === #
latent_vars = list(
# === Non linear model === #
latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
# === S.data model === #
latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
add_on = c('zeta(b = b, ...) +',
'sum(h$eval(b = b, ..., i = c(1,2)))' )),
latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
add_on = c('zeta(alpha = alpha, ...) +',
'alpha*h$eval(alpha = alpha,..., i = 3)'))
),
# === Paramètre de regression === #
regression.parameter = list(
regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
step = 0.05, lambda = 1/sqrt(N),
normalized.grad = T,
zeta.der.B, N, zeta.B,
Z$alpha, Z$phi1, Z$phi2, Z$phi3,Z$b) )
)
)
# --- Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))
#===============================================#
load.SAEM(model)
U <- U_full
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#
init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha),
sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )
SAEM.options <- list(niter = 300, sim.iter = 5, burnin = 190,
adptative.sd = 0.6)
saem
p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 32min 38sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9034
|
89.9207
|
5.0293
|
0.0039
|
35.3106
|
0.6393
|
28.7713
|
0.8068
|
|
Rrmse
|
0.0008
|
0.0002
|
0.0004
|
0.0043
|
0.0194
|
0.0169
|
0.0798
|
0.0410
|
0.6135
|
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 36min 24sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9045
|
89.9446
|
5.0482
|
0.0039
|
35.5880
|
0.6418
|
33.5285
|
0.6676
|
|
Rrmse
|
0.0023
|
0.0015
|
0.0001
|
0.0081
|
0.0152
|
0.0092
|
0.0762
|
0.1176
|
0.3352
|
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 42min 04sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9034
|
89.9501
|
5.0219
|
0.0039
|
35.0544
|
0.6071
|
36.9021
|
2.2173
|
|
Rrmse
|
0.0012
|
0.0002
|
0.0001
|
0.0028
|
0.0199
|
0.0240
|
0.1262
|
0.2301
|
3.4346
|
p <- 500
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 37min 12sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9028
|
89.9222
|
5.0191
|
0.0038
|
34.9631
|
0.6428
|
25.2145
|
1.1898
|
|
Rrmse
|
0.0022
|
0.0004
|
0.0004
|
0.0023
|
0.0431
|
0.0266
|
0.0748
|
0.1595
|
1.3797
|
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 34min 15sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9030
|
89.9085
|
5.0218
|
0.0039
|
35.1083
|
0.6076
|
22.1310
|
0.5723
|
|
Rrmse
|
0.0006
|
0.0003
|
0.0005
|
0.0028
|
0.0100
|
0.0225
|
0.1255
|
0.2623
|
0.1447
|